function H = cr_matrix(dofs,V,T,TE,ET,desc,asce,cr,d,r)
% function H = cr_matrix(dofs,V,T,TE,ET,desc,asce,cr,d,r)
% get the global smoothness conditions
% Assuming that the degree d is vector!!

innerEdges = find(ET(:,2)>0);
nE = length(innerEdges);
I_H1=zeros(nE*(4*max(d)+1),1); J_H1=I_H1; S_H1=I_H1;
I_H2=zeros(20*nE*(2*max(d)+1),1); J_H2=I_H2; S_H2=I_H2;
M_bary = ones(3,3);
b_bary = ones(3,1);
row = 1; pos1 = 1; pos2 = 1;
for k = 1:nE
    eg = innerEdges(k);
    me = ET(eg,1); nei = ET(eg,2);
    if d(nei) > d(me) % always keep the degree of me is the bigger;
        tmp = me;
        me = nei;
        nei = tmp;
    end
    degree = d(me);
    % begin to treat the matrix H1 which get from triangle me
    % evaluate barycentric coordinate of v4:
  if r == 1
    v4 = find(TE(nei,:)==eg);  
    M_bary(2,1) = V(T(me,1),1); M_bary(2,2) = V(T(me,2),1); M_bary(2,3) = V(T(me,3),1);
    M_bary(3,1) = V(T(me,1),2); M_bary(3,2) = V(T(me,2),2); M_bary(3,3) = V(T(me,3),2);
    b_bary(2) = V(T(nei,v4),1);  b_bary(3) = V(T(nei,v4),2);
    mu = M_bary\b_bary;
  end
    eg_in_me = find(TE(me,:)==eg);
    if eg_in_me == 2 % eg in me is 2 ,using anti-order because the eg is v3v1
        eg_in_me = -2;
    end
    % for c^0 condition:
    line_me = cr_indices(0,degree,eg_in_me,cr);
    I_H1(pos1:pos1+degree) = row:row+degree;
    J_H1(pos1:pos1+degree) = dofs(me,1) - 1 + line_me;
    S_H1(pos1:pos1+degree) = ones(degree+1,1);
    pos1 = pos1 + degree + 1;
    % for c^1 condition, we have an alternative way to do this! 
    if r == 1 
        line_me = cr_indices(0,degree-1,eg_in_me,cr);
        mat_desc = desc_mat(degree,mu(1),mu(2),mu(3),desc); 
        mat_desc = mat_desc(line_me,:);
        [i,j,s] = find(mat_desc);
        L = length(i);
        I_H1(pos1:pos1+L-1) = row + degree + i;
        J_H1(pos1:pos1+L-1) = dofs(me,1) - 1 + j;
        S_H1(pos1:pos1+L-1) = s;
        pos1 = pos1 + L;    
    end    
    % treate the matrix H2:
    eg_in_nei = find(TE(nei,:)==eg);
    if eg_in_nei ~= 2 % eg in me is not 2 ,using anti-order because the eg is v2v1 or v3v2
        eg_in_nei = -1*eg_in_nei;
    end
    line_nei = cr_indices(0,degree,eg_in_nei,cr);
    if r==1 % add the c^1 condition line's dof
        line_nei = [line_nei;cr_indices(1,degree,eg_in_nei,cr)];
    end
    %then get the elevate matrix to the same degree as me
    if degree > d(nei) % only the degrees are different, we need to do this 
        mat_asce = asce_mat(degree,asce);
        mat_asce = mat_asce(line_nei,:); % only get the needed dofs
        for di = (degree-1):-1:d(nei)+1 % may const some time
            mat_asce = mat_asce * asce_mat(di,asce);
        end
        [i,j,s] = find(mat_asce);
        L = length(i);
        I_H2(pos2:pos2+L-1) = row - 1 + i;
        J_H2(pos2:pos2+L-1) = dofs(nei,1) - 1 + j;
        S_H2(pos2:pos2+L-1) = s;
        pos2 = pos2 + L;
    elseif degree == d(nei) % the identity matrix
        L = length(line_nei);
        I_H2(pos2:pos2+L-1) = (row:row+L-1)';
        J_H2(pos2:pos2+L-1) = dofs(nei,1) - 1 + line_nei;
        S_H2(pos2:pos2+L-1) = ones(L,1);
        pos2 = pos2 + L;
    end
    row = row + length(line_nei);
end
H = sparse(I_H1(1:pos1-1),J_H1(1:pos1-1),S_H1(1:pos1-1),row-1,dofs(end,2)) - sparse(I_H2(1:pos2-1),J_H2(1:pos2-1),S_H2(1:pos2-1),row-1,dofs(end,2));